// StepProbInterface.C

#include <iostream>
#include <fstream>
#include <string>

#include "Epetra_CrsMatrix.h"
#include "LOCA_Parameter_Vector.H"

#include "StepProbInterface.H"
#include "StepWaveProb.H"
#include "Global.H"

StepProbInterface::StepProbInterface(StepWaveProb &prob, char *contparam)
  : contParam(contparam),
    prevEval(EVALINIT),
    problem(prob)
{ 
}

StepProbInterface::~StepProbInterface()
{ }

void StepProbInterface::printSolution(const::Epetra_Vector &x, double paramVal)
{
  cerr << contParam << " = " << paramVal << endl;

  int m = problem.getM();
  int n = problem.getN();
  int dim = m * n;
  double p0 = problem.getP0();
  double q = problem.getQ();
  double p1 = problem.getP1();
  double p2 = problem.getP2();
  double a1 = problem.getA1();
  double a2 = problem.getA2();
  double a3 = problem.getA3();

  double avg = 0.0;
  for (int i=0; i<m; ++i)
    avg += x[i];
  avg /= m;

  double d = 0.0;
  for (int i=0; i<m; ++i)
    d += (avg - x[i])*(avg - x[i]);

  cerr << "  sigma = " << sqrt(d/n) << endl;


  char filename[32];
  sprintf(filename, "sol.%s%.4f", contParam, paramVal);
  ofstream fout(filename, ios::trunc);
  fout.precision(15);
  fout << "step" << endl
       << m << endl
       << n << endl
       << p0 << endl
       << q << endl
       << p1 << endl
       << p2 << endl
       << a1 << endl
       << a2 << endl
       << a3 << endl << endl;
  for (int i=0; i<dim; ++i)
    fout << x[i] << endl;
  fout.close();
}

void StepProbInterface::printEigenvector(const::Epetra_Vector &x, double eval, int which)
{
  cerr << "    eval " << which << " = " << eval << endl;

  if (which != 1) return;

  if (prevEval == EVALINIT)
    prevEval = eval;
  else if (prevEval*eval < 0)
  {
    cerr << "Detected sign change in second eigenvalue; writing " << endl
         << "eigenvector to file and stopping continuation." << endl;
    char filename[32];
    sprintf(filename, "evec%d.%4f", which, eval);
    //sprintf(filename, "evec%d", which);
    ofstream fout(filename, ios::trunc);
    fout.precision(15);
    int dim = problem.getM() * problem.getN();
    for (int i=0; i<dim; ++i)
      fout << x[i] << endl;
    fout.close();
   // exit(1);
  }
  else
  prevEval = eval;
}

bool StepProbInterface::computeF(const Epetra_Vector& x, Epetra_Vector& F, FillType ft)
{
  return problem.computeF(x, F);
}

bool StepProbInterface::computeJacobian(const Epetra_Vector& x)
{
  return problem.computeJacobian(x);
}

void StepProbInterface::setParameters(const LOCA::ParameterVector& params)
{
  char label[2];
  double value;
  for (int i=0; i<params.length(); ++i)
  {
    sprintf(label, "%s", params.getLabel(i).c_str());
    value = params.getValue(i);
    if (strcmp(label, "Q") == 0)
      problem.setQ(value);
    else if (strcmp(label, "P1") == 0)
      problem.setP1(value);
    else if (strcmp(label, "P2") == 0)
      problem.setP2(value);
    else if (strcmp(label, "A1") == 0)
      problem.setA1(value);
    else if (strcmp(label, "A2") == 0)
      problem.setA2(value);
    else if (strcmp(label, "A3") == 0)
      problem.setA3(value);
  }
}

